1. Field of the Invention
The invention relates to signal processing systems and methods in general and particularly to a signal processing method that employs recursive interpolation and sparse-matrix methods implementation with thresholding in order to speed up the search for the parameters of a linear frequency modulated signal.
2. Description of the Prior Art
The linear frequency modulated (FM) signal or “chirp” is used extensively in signal processing to model man-made or natural signals. When the parameters (starting frequency and slope or chirp rate) are unknown, it is necessary to search a two-dimensional parameter space. Since there are on the order N2 useful pairs of starting frequencies and slopes to test, a brute-force search would require a computational load of N3 for a length-N input data record. The state-of-the-art method using the Fast Fourier Transform (“FFT”) achieves a computational load of N2 log N.
Consider the linear frequency-modulated (FM) chirp signal defined bysi(a,γ,f,θ)=a exp{jti2γ+j2πti+jθ},1≦i≦N,  (1)where N and i are integers, a is an amplitude, θ is a phase angle parameter, the slope parameter γ exists in the rangeNπ/2≦γ≦Nπ/2the frequency parameter f exists in the range0≦f≦N/2,and the time samples areti=(i−1)/N. 
When γ,f are unknown, the maximum likelihood (ML) solution searches over the two dimensional space (γ,f). Let γ be selected from the set
      _    ⁢                  ⁢    Γ    =                                  (                      -                          N              2                                )                ⁢        π            ,                        (                                    -                              N                2                                      +            δ                    )                ⁢        π            ,              …        ⁢                                  ⁢                  (                      N            2                    )                ⁢        π                where δ is the slope search resolution parameter. If δ=2, there are N/2+1 possible slope values. Since we have found that δ=2 provides adequate coverage, we will assume from here on that δ=2.
Let f be selected from the set
  F  =                          0        ,                  1          2                ,        1        ,                  3          2                ,                  …          ⁢                                          ⁢                      N            2                                      .  
Note that the frequency search resolution is 2 times finer than an FFT requiring that a size-2N transform be used. There are (N+1)(N/2+1) bins to search for this sampling of the (γ,f) plane. The search surface function is given by:
                              F          ⁡                      (                                                            x                  1                                ⁢                                                                  ⁢                …                ⁢                                                                  ⁢                                  x                  N                                            ,              f              ,              γ                        )                          =                              ∑                          i              =              1                        N                    ⁢                                    x              i                        ⁢            exp            ⁢                                          {                                  -                                      j                    ⁡                                          (                                                                                                    t                            i                            2                                                    ⁢                          γ                                                +                                                  2                          ⁢                          π                          ⁢                                                                                                          ⁢                                                      t                            i                                                    ⁢                          f                                                                    )                                                                      }                            .                                                          (        2        )            
The brute-force calculation of the search surface |F(x1 . . . xN, f, γ)| is implemented in MATLAB®, available from The MathWorks, Inc., 3 Apple Hill Drive, Natick, Mass. 01760-2098. The MATLAB code is presented in TABLE I below and shows the elapsed time for the calculation.
TABLE I%  set data sizeN=256;%  choose injected signal parameters. gam is slope, frq is  frequencygam=−N/8*pi; frq=N/4;%  set search bins for slope and frequencygams=[−N/2 : 2 : N/2]*pi;frqs=[0: 1/2 : N/2]′;%  create signal and noiseSIG_LEVEL=1.0;t=[0:N−1]′/N;s=exp(1i*t.{circumflex over ( )}2*gam+1i*2*pi*t*frq);x=SIG_LEVEL * sqrt(128/N) * real(s)+randn(N,1);%  brute-force search%  fill array with zerosF=zeros(N+1,N/2+1);tic%  loop from 1 to N/2 +1 where j is the loop variable for  slope  for j=1:N/2+1,%  loop from 1 to N+1 where i is the loop variable for  frequency   for i=1:N+1,%  calculate periodic portion of the search surface    e=exp(1i*t.{circumflex over ( )}2*gams(j)+1i*2*pi*t*frqs(i));%  multiply periodic portion by coefficients   F(i,j) = e′*x;  end; end;tocelapsed_time = 5.9481imagesc(gams,frqs,abs(F)); axis xyxlabel (‘slope (gamma)’);ylabel (‘frequency (f)’);
FIG. 1 is a diagram that illustrates a typical search surface.
Since correlating the input data at each test frequency and slope requires N operations, the process is of order N(N+1)(N/2+1)≅N3/2. This is prohibitive if N is large.
The maximum likelihood search is typically done using the fast Fourier transform (FFT) to greatly speed up the search over f. For a given choice of γ, we form the de-chirped signal
            z      i        =                  x        i            ⁢      exp      ⁢              {                              -            j                    ⁢                                          ⁢                      t            i            2                    ⁢          γ                }              ,where xi are input data samples containing signal si plus additive independent Gaussian noise. Because the de-chirped signal needs to be zero-padded to length 2N, the total computation load for the FFT is 2 Nlog2(2N), which must be computed N/2+1 times for a total computational load of (N/2+1)2 Nlog2(2N).
A reduction is possible when taking advantage of the zero-padded FFT, bringing the computational load to (N/2+1)2Nlog2(N). The MATLAB code is presented in TABLE II below.
TABLE IIE=zeros(N,N/2+1);for j=1:N/2+1, E(:,j)=exp(−1i*t.{circumflex over ( )}2*gams(j));end;ticX=repmat(x,1,N/2+1);F=fft(X.*E,2*N);tocelapsed_time = 0.0115imagesc(gams,frqs,abs(F)); axis xyxlabel (‘slope (gamma)’);ylabel (‘frequency (f)’);
The MATLAB code given in TABLE II produces exactly the same surface as the brute-force method. Our observed speedup factor of about 550 is mainly due to the vectorization that we used in the FFT approach as opposed to the inefficient looping that we used in the brute-force algorithm. Despite the improvements, the order N2 log(N) computational load still too high for many applications.
Known in the prior art is Gruber et al., U.S. Pat. No. 5,249,150, which is said to disclose a circuit and method for iteratively estimating Fourier coefficients of one or more Fourier components of a measuring signal that utilizes a filter to determine the Fourier coefficients for sampling time k−1. The coefficients are utilized to obtain an estimate of the measuring signal at sampling time k. A subtractor circuit subtracts an actual sample of the measuring signal at sampling time k from the estimate of the measuring signal at sampling time k to obtain a difference signal. The difference signal is inputted back into the filter to determine an estimate of the Fourier coefficients at sampling time k.
Known in the prior art is Fukuda et al., U.S. Pat. No. 5,960,373, which is said to disclose a frequency analyzing method for analyzing frequency components of an original signal. The method has: a spectrum detecting step of detecting, from the original signal, energy levels of components of a predetermined number of orthogonal function waves which have waveforms each having same start position and end position in a predetermined time window and in which the number of occurrences of periods in the predetermined time window or frequencies are different from each other; and an orthogonal function wave changing step of changing at least one of the start position and the end position within the predetermined time window after completion of the spectrum detecting step, wherein the spectrum detecting step and the orthogonal function wave changing step are alternately repeated. According to Fukuda, it is possible to provide frequency analyzing method and apparatus which can contribute to estimate each correct fundamental frequency from a complex distorted wave signal such as a musical signal or the like by a relatively simple construction and to provide complex sound separating method and apparatus using the frequency analyzing method or apparatus.
Also known in the prior art is Scheppach, U.S. Pat. No. 6,484,112, which is said to disclose a method for estimating the frequency of a time signal by means of a discrete Fourier transformation (DFT) and interpolation between support points of the DFT spectrum. According to the method, the Hamming window is used in a known manner for filtering and the interpolation is carried out according to the mathematically calculated solution. The method makes use of the fact that a function of the third degree which can serve as a correction distance has a clear solution within the definition range of the amount ratio between a secondary maximum and a primary maximum of the DFT spectrum. The correction distance is calculated according to this solution on the basis of the ratio of the secondary maximum and the primary maximum.
Also known in the prior art is Nelson, U.S. Pat. No. 6,577,968, which is said to disclose a method of estimating the frequency of a signal by first receiving a signal, forming a row vector, segmenting the row vector; converting the row vector to a first matrix, multiplying the first matrix by a weight, and performing a Fourier transform on the result of the last step. These same steps are repeated on a delayed version of the signal. The next steps are calculating a complex conjugate for each result of the last step, forming a cross-spectrum matrix, selecting a magnitude in the cross-spectrum matrix that is above a threshold; and setting an angular frequency of the signal to either the phase of the selected magnitude, the phase of the mean of the complex numbers in the row in which appears the selected magnitude, or the selected magnitude. The frequency of the signal is then set to the estimated angular frequency divided by the product of 2π and the signal delay period. Additional rows are processed if need be.
Also known in the prior art is Hillerstrom, U.S. Pat. No. 6,598,005, which is said to disclose a method for measuring a frequency f=ω/2π of a sinusoidal signal V(t)=A sin(ωt+α) with an essentially constant amplitude A. Measured values are registered in the form of the instantaneous level of the signal V(t), V(t−h) and V(t−2 h) at three points of time separated by a predetermined measuring period h, which corresponds to a measurement frequency fs=1/h. Then the measured values are applied to a calculating device, which recursively estimates the quantity y(t) given by: y(t)=y(t−h)+ηe(t)s(V(t−h)), where e(t)=[V(t)+V(T−2 h)]/2−y(t−h) and ηs(V(t−h))V(t−h)≧0. The choice of η and s(.) controls the convergence speed of y(t)→cos(ωh) and the noise sensitivity of the measurement system. Finally the frequency is calculated or looked up in a table based on the expression f(t)=(fs/2π)arc cos y(t).
Also known in the prior art is Michel, U.S. Pat. No. 7,076,380, which is said to disclose a method for analyzing the frequency in real-time of a non-stationary signal and corresponding analysis circuit. According to the invention, the signal is sampled, digitized, and broken up into sub-bands. In each sub-band, the signal is modeled by an auto-regressive filter the transfer function of which is of the form 1/A(z). By an adaptive method that is recursive in time and in order, all the polynomials A(z) that have a degree between 1 and a maximum value are calculated. The order of the model is estimated and the polynomial that has this order is retained. The roots of this polynomial are calculated and the components are monitored. The frequency and the amplitude of the sinusoidal components of the signal are thus obtained. Applications include aeronautics, electromagnetic, mechanics, seismic prospecting, zoology, and others.
Also known in the prior art is Shim et al., U.S. Pat. No. 7,103,491, which is said to disclose a method of estimating parameters of time series data using a Fourier transform. The parameter estimation method can directly estimate modes and parameters, indicating the characteristics of a dynamic system, on the basis of the Fourier transform of the time series data. The method is advantageous in that it can estimate parameters of the system through a single Fourier transform without repeatedly calculating a Fourier transform, thus shortening the time required to estimate the parameters.
Also known in the prior art is Rao, U.S. Pat. No. 7,124,042, which is said to disclose a system and method for estimating parameters of multiple tones in an input signal. The method includes receiving samples of the input signal, generating a frequency transform (FT) of the samples, identifying multiple amplitude peaks in the FT corresponding to the tones, and determining parameter estimates characterizing each of the multiple tones based on the peaks. For each tone, the effects of the other tones are removed from the FT of the peak of the tone using the parameter estimates of the other tones to generate modified FT data for the tone. Single tone estimation is applied to the modified FT data to generating refined parameter estimates of the tone, which is used to update the parameter estimates of the tone. After refining the estimates for each tone, the entire process may be repeated one or more times using successive refined estimates to generate final estimates for the parameters.
A literature search revealed a number of algorithms related to this problem that include the following. There are algorithms that attain fast performance by searching only over 1-dimensional space of γ (see J. Cao, N. Zhang, and L. Song, “A fast algorithm for the chirp rate estimation,” in Proceedings of the 4th IEEE International Symposium on Electronic Design, Test and Applications (DELTA 2008), pp. 45-48, January 2008; P. M. Djuric and S. M. Kay, “Parameter estimation of chirp signals,” IEEE Trans on Acoustics, Speech and Signal Processing, vol. 38, pp. 2118-2126, December 1990; and W. Li, M. Dan, X. Wang, D. Li, and G. Wang, “Fast estimation method and performance analysis of frequency modulation rate via rat,” in Proceedings of the 2008 International Conference on Information and Automation (ICIA), (Zhangjiajie, China), pp. 144-147, June 2008). These algorithms confuse signals with different starting frequency and similar slope and typically require a higher signal-to-noise ratio (SNR) to work properly.
There are algorithms that integrate in the time-frequency transform (TF) plane (see S. Barbarossa, “Analysis of multicomponent lfm signals by a combined wigner-hough transform,” IEEE Trans on Signal Processing, vol. 43, pp. 1511-1515, June 1995; S. Krishnan and R. Rangayyan, “Detection of chirp and other components in the time-frequency plane using the hough and radon transforms,” in Proceedings of the 1997 IEEE Pacific Rim Conference on Communications, Computers and Signal Processing, vol. 1, (Victoria, BC, Canada), pp. 138-141, August 1997; and J. Wood and D. Barry, “Radon transformation of time-frequency distributions for analysis of multicomponent signals,” IEEE Trans on Signal Processing, vol. 42, pp. 3166-3177, November 1994). The computation of the TF distribution adds significant processing load. In addition, because of cross-terms, these methods can have difficulties with multi-component signals.
There are algorithms that are based on the fractional Fourier transform (FrFT) See H. M. Ozaktas, O. Ankan, M. A. Kutay, and G. Bozdaki, “Digital computation of the fractional fourier transform,” IEEE TRANSACTIONS ON SIGNAL PROCESSING, vol. 44, September 1996. No significant processing improvement can be obtained because of the necessity of calculating a FrFT for each assumed slope.
A set of similar algorithms called various names such as the discrete chirp Fourier transform (DCFT) (See C. A. Aceros-Moreno and D. Rodriguez, “Fast discrete chirp fourier transforms for radar signal detection systems using cluster computer implementations,” in Proceedings of the 48th Midwest Symposium on Circuits and Systems, vol. 2, pp. 1047-1050, August 2005), discrete quadratic phase transform (DQPT) (See M. Z. Ikram, K. Abed-Meraim, and Y. Hua, “Fast discrete quadratic phase transform for estimating the parameters of chirp signals,” Digital Signal Processing, vol. 7, pp. 127-135, April 1997), and polynomial time-frequency transform (PTFT) (See G. Bi and Y. Wei, “Split-radix algorithms for arbitrary order of polynomial time frequency transforms,” IEEE Trans on Signal Processing, vol. 55, pp. 134-141, January 2007; and Y. Ju and G. Bi, “Generalized fast algorithms for the polynomial time-frequency transform,” IEEE Trans on Signal Processing, vol. 55, pp. 4907-4915, October 2007). The DCFT is defined as:
            X      ⁡              [                  k          ,          l                ]              =                  ∑                  i          =          0                          N          -          1                    ⁢                        x          i                ⁢                  W          N                                    li              2                        +            ki                                ,where WN=ej2π/N, l and k are integers in the range 0 to N−1. Since useful chirp signals use fractional values of l between −½ and ½, this transform that uses integer l is of little value for our purposes. The authors claim an impressive computational reduction, but this requires the calculation of mostly unrealistic values of slope. The DQPT and PTFT are similarly defined and have the same problems.
In short, one is left with the FFT-based algorithm described above that has N2 log2(N) computational load as the only realistic solution available up to now. There is a need for more efficient and less time consuming methods for analyzing linear frequency-modulated signals.